Assignment 3 - The free position facility location problem

The Problem

In this assignment we will look at the Facility Location Problem. The problem can be described in general terms as follows:

We must decide the placement of a facility that distributes a certain good to a group of consumers that need it.

The placement must to be chosen in order to minimize the total compound distance from the facility and the customers.

The following assumptions has to be taken into account:

1. The possible location for the facility is unknown, that is the problem is to find the right spot to build it.
2. The facility building costs are fixed and independent from the position of the building site.

Notice that in this scenario there is one possible decision to make:

• where to build the facility, that is find the position (\(\chi\), \(\upsilon\)) that minimises the compound distance of the facility with respect to all custumers.

Distance function

Given the position of the facility \(f= (\chi, \upsilon)\) and of a consumer \(p_i=(x_i, y_i)\) use the following formula to calculate the distance between them.

\[ d(f,p_i) = log((\chi-x_i)^2+1) + log((\upsilon-y_i)^2+1) \] #Carico le coordinate su una variabile:

library(ggplot2)
coordinate<-read.csv("/Users/lorenzofamiglini/Downloads/consumer_locations.csv")

1. Formulate the objective function to minimize for the described problem. La mia funzione obiettivo non è altro che la minimizzazione della funzione del problema:

\[ min \ f(\chi, \upsilon) \\ \text{ dove }f(\chi,\upsilon)=\sum_{i=1}^{n} log((\chi-x_i)^2+1) + log((\upsilon-y_i)^2+1) \] Le derivate parziali sono: \[ ∂/∂\chi = \sum_{i=1}^{n}(2(\chi-x_i)/(\chi-x_i)^2 + 1) \\ ∂/∂\upsilon = \sum_{i=1}^{n}(2(\upsilon-y_i)/(\upsilon-y_i)^2 + 1) \] Rappresentiamo la funzione in 3D:

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

plottare <- function(x, y){
  distanza <- 0 
  for (i in 1:nrow(coordinate)){
  distanza <- distanza + log((x-coordinate$x[i])^2+1) + log((y - coordinate$y[i])^2 + 1)
  }
  return(distanza)
}

x <- seq(0, 1000)
y <- x
z = outer(x, y, plottare)
plt <- plot_ly(z = ~ z, colors = c("blue","grey")) %>% add_surface(
  contours = list(
    z = list(
      show = TRUE,
      usecolormap = TRUE,
      highlightcolor = "#ff0001",
      project = list(z = TRUE)
      )
    )
  ) %>%
  layout(
    scene = list(
      camera = list(
        eye = list(x = 1.87, y = 0.88, z = - 0.64)
        )
      )
  )
plt

Creiamo una funzione che calcola la funzione obiettivo in base a dei punti (x,y):

# f<- expression(log((x-xi)^2+1) + log((y - yi)^2 + 1)) #l'espressione generale
func_loss <- function(punto){
  x <- punto[1] #primo elemento del vettore
  y <- punto[2] #secondo elemento del vettore 
  sum(log((x-coordinate[1])^2+1) + log((y - coordinate[2])^2 + 1))
}
func_loss(c(40,40)) #ad esempio mettiamo le coordinate nel punto(40,40) 

## [1] 2273.076

2. Express in analytical form the gradient for the objective to minimize.

#Calcoliamo le derivate parziali rispetto ad un punto (x,y) creando una funzione che implementa tale operazione:

pard <- function(punto, coordinate){ #partial derivates (pard)
  m <- as.matrix(coordinate)
  p1 <- punto[1]
  p2 <- punto[2]
  der1 <- 0 
  der2 <- 0 
  for (i in 1:nrow(m)) {
    der1 <- der1 + 2*(p1 - m[i,1])/((p1-m[i,1])**2 + 1) #derivata parziale rispetto a chi
    der2 <- der2 + 2*(p2 - m[i,2])/((p2-m[i,2])**2 + 1) #derivata parziale rispetto a upsilon
  }
  return(c(der1,der2))
  
}
pard(c(40,40), coordinate)

##          x          y 
## -1.7245915  0.7016332

Questa funzione verrà implementata nei successivi algoritmi.

3. Implement the Gradiend Descent method and solve the problem with it.

gdm<- function(coordinate, punto, parameter, n_iter, func_loss) {
  it <- 1
  l <- c()
  v <- c()
  p <- c()
  l[it] <- it
  v[it] <- func_loss(punto)
  p[it] <- parameter
  step <- pard(punto, coordinate)*parameter
  while (any((abs(step)) > c(0.00001, 0.00001)) && it <=n_iter)
    {
    step <- pard(punto, coordinate)*parameter #derivata parziale per il learning rate
    punto <- punto - step #nuove coordinate del punto
    parameter <- parameter*0.95
    it <- it +1
    l[it] <- it
    v[it] <- func_loss(punto)
    p[it] <- parameter
  }
  list(minimo = v, Coordinate = punto,
       iterazioni = l, learning_rate = p)
}

La funzione è stata implementata secondo l’algoritmo del gradient descent method, all’interno di essa abbiamo un ciclo while che si ferma in due occasioni: se il numero di iterazioni ha raggiunto il massimo oppure se lo step size è molto vicino allo zero. E’ stato utilizzato un parametro che diminuisce all’aumentare delle iterazioni, senza lasciarlo costante. Infatti, ogni volta che finisce l’iterazione il learning rate è moltiplicato per 0.9, in modo tale da creare “passi” sempre più piccoli. Le coordinate nuove vengono calcolate sottraendo il punto vecchio con lo step size.

Mettiamo a confronto vari parametri iniziali, prendendo un punto iniziale casuale (40,40):

result1 <- gdm(coordinate, c(40,40), 1, 100, func_loss)
result2 <- gdm(coordinate, c(40,40), 2, 100, func_loss)
result3 <- gdm(coordinate, c(40,40), 3, 100, func_loss) 
result1 <- data.frame(distanza = result1$minimo, iterazioni = result1$iterazioni, learning_rate =result1$learning_rate)
result2 <- data.frame(distanza = result2$minimo, iterazioni = result2$iterazioni, learning_rate =result2$learning_rate)
result3 <- data.frame(distanza = result3$minimo, iterazioni = result3$iterazioni, learning_rate =result3$learning_rate)

result1

result2

result3 

ggplot() +
  geom_line(data = result1, aes(x = iterazioni, y=distanza,  color = "red")) + 
  geom_line(data = result2, aes(x = iterazioni, y=distanza, color = "blue")) + 
  geom_line(data = result3, aes(x = iterazioni, y=distanza, color = "green")) + 
   labs(title="Andamento dell'algoritmo per ogni iterazione", 
        y = "distanza", 
        x = "numero di iterazioni") +
    scale_color_manual("Learning rate \n di partenza", values = c("red", "blue", "green"), labels = c("1","2","3"))

Il seguente grafico mostra che partire da un parametro più grande (2) permette di raggiungere un punto di minimo migliore rispetto ad inizializzare l’algoritmo con parametri più piccoli. Questa considerazione ha validità partendo da un punto di coordinate (40,40). La distanza ottimale individuato da tale funzione è: 2269.407 (molto probabilmente sarà un minimo locale, alla fine dell’analisi andiamo a cercare il punto ottimale di partenza).

4. Solve the problem with a package provided by R (for instance, using the function optimr within the package optimx). Note that it is not required to use the gradient descent algorithm to solve the problem, other algorithms can be used as well.

library(optimx)
ss <- optimr(c(40,40), func_loss, method="CG") #cognitive gradient
ss

## $par
## [1] 42.23663 39.10584
## 
## $value
## [1] 2270.518
## 
## $counts
## function gradient 
##       21       11 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

proptimr(ss) # compact output of result

## Result  ss   proposes optimum function value = 2270.518  at parameters
## [1] 42.23663 39.10584
## After  21  fn evals, and  11  gr evals
## Termination code is  0 :
## ---------------------------------------

Utilizzando la libreria optmix, si raggiungo dei risultati pressochè simili all’algoritmo sviluppato in modo analitico, anzi sembrerebbe che la funzione precedente sia leggermente più precisa.

5. Implement the Stochastic gradient descent algorithm with mini-batches and use it to solve the problem.

sgdm<- function(coordinate, punto, parameter, n_iter, func_loss) {
  it <- 1
  v <- c()
  l <- c()
  p <- c()
  step <- pard(punto, coordinate)*parameter
  v[it] <- func_loss(punto)
  l[it] <- it
  p[it] <- parameter
  while (any((abs(step)) > c(0.00001, 0.00001)) && it <= n_iter)
    {
    set.seed(123)
    mini_batch <- coordinate[sample(nrow(coordinate), 15),] #prendiamo il 15% delle osservazioni dal dataset
    step <- pard(punto, mini_batch)*parameter #derivata parziale per il learning rate
    punto <- punto - step #nuove coordinate del punto
    parameter <- parameter*0.95
    it <- it +1 
    l[it] <- it
    v[it] <- func_loss(punto)
    p[it] <- parameter
  }
  list(minimo = v, Coordinate = punto,
       iterazioni = l, learning_rate = p)
}

Per implementare tale algoritmo, la logica alla base è identica al classico GDM, cambia semplicemente il fatto che si utilizzano dei campioni (si è scelto il 15% del totale delle osservazioni) per poter calcolare i vari punti. In questo caso il metodo non stocastico è più accurato, anche perchè possiamo sfruttare tutto il campione (avendo un’ampiezza molto piccola). Ma nel caso in cui ci troviamo di fronte ad una grande mole di dati è molto consigliato usare il GDM stocastico.

#Testiamo il parametro ottimale di partenza e vediamo come è l’andamento fino a raggiungere il punto di ottimo (partendo da un punto, si ribadisce ancora una volta che non è detto che sia globale, anzi è più probabile che sia locale). Il parametro iniziale cambia per ogni risultato. Questo ci fa capire quale è il parametro iniziale migliore per trovare il punto di minimo. Facciamo un esempio partendo da un punto iniziale (40,40):

result1 <- sgdm(coordinate, c(40,40), 1, 100, func_loss)
result2 <- sgdm(coordinate, c(40,40), 2, 100, func_loss)
result3 <- sgdm(coordinate, c(40,40), 3, 100, func_loss) 
result4 <- sgdm(coordinate, c(40,40), 4, 100, func_loss) 
result5 <- sgdm(coordinate, c(40,40), 5, 100, func_loss) 
result6 <- sgdm(coordinate, c(40,40), 6, 100, func_loss)
result7 <- sgdm(coordinate, c(40,40), 7, 100, func_loss)
result8 <- sgdm(coordinate, c(40,40), 8, 100, func_loss)


result1 <- data.frame(distanza = result1$minimo, iterazioni = result1$iterazioni, learning_rate =result1$learning_rate)
result2 <- data.frame(distanza = result2$minimo, iterazioni = result2$iterazioni, learning_rate =result2$learning_rate)
result3 <- data.frame(distanza = result3$minimo, iterazioni = result3$iterazioni, learning_rate =result3$learning_rate)
result4 <- data.frame(distanza = result4$minimo, iterazioni = result4$iterazioni, learning_rate =result4$learning_rate)
result5 <- data.frame(distanza = result5$minimo, iterazioni = result5$iterazioni, learning_rate =result5$learning_rate)
result6 <- data.frame(distanza = result6$minimo, iterazioni = result6$iterazioni, learning_rate =result6$learning_rate)
result7 <- data.frame(distanza = result7$minimo, iterazioni = result7$iterazioni, learning_rate =result7$learning_rate)
result8 <- data.frame(distanza = result8$minimo, iterazioni = result8$iterazioni, learning_rate =result8$learning_rate)

ggplot() +
  geom_line(data = result1, aes(x = iterazioni, y=distanza, color = "red")) + 
  geom_line(data = result2, aes(x = iterazioni, y=distanza, color = "blue")) + 
  geom_line(data = result3, aes(x = iterazioni, y=distanza,  color = "green")) + 
  geom_line(data = result4, aes(x = iterazioni, y=distanza, color = "orange")) + 
  geom_line(data = result5, aes(x = iterazioni, y=distanza,  color = "purple")) + 
  geom_line(data = result6, aes(x = iterazioni, y=distanza, color = "black")) + 
  geom_line(data = result7, aes(x = iterazioni, y=distanza, color = "tomato")) + 
  geom_line(data = result8, aes(x = iterazioni, y=distanza, color = "grey")) +
    labs(title="Andamento dell'algoritmo per ogni iterazione", 
        y = "distanza", 
        x = "numero di iterazioni") +
  scale_color_manual("Learning rate \n di partenza", values = c("red", "blue", "green", "orange", "purple", "black","tomato","grey"), labels = c("1","2","3", "4","5","6","7","8"))

result2

Il grafico mostra che partendo da un learning rate molto basso (2 circa), riusciamo a raggiungere un punto di minimo migliore (con una distanza pari a:2270.518).Ovviamente sono stati utilizzati diversi learning rate di partenza per vedere se l’algoritmo saltava il punto di minimo per cascare in uno errato. Infatti, si nota che se il “passo” è troppo lungo saltando il candidato punto di minimo ottimale. NB: I risultati cambieranno in base al punto di partenza, questo è un grafico puramente a scopo “illustrativo”

Il successivo grafico mostra l’andamento del learning rate all’aumentare delle iterazioni: Il parametro decresce per ogni iterazione. Infatti, esso è moltiplicato per un valore vicino ad uno (0.9) in modo tale da poter “rallentare” l’andamento.

ggplot() +
  geom_line(data = result1, aes(x = iterazioni, y=learning_rate), color = "red") +
  geom_line(data = result2, aes(x = iterazioni, y=learning_rate), color = "blue") + 
  geom_line(data = result3, aes(x = iterazioni, y=learning_rate), color = "green") + 
  geom_line(data = result4, aes(x = iterazioni, y=learning_rate), color = "orange") + 
  geom_line(data = result5, aes(x = iterazioni, y=learning_rate), color = "purple") + 
  geom_line(data = result6, aes(x = iterazioni, y=learning_rate), color = "black") + 
  geom_line(data = result7, aes(x = iterazioni, y=learning_rate), color = "tomato") + 
  geom_line(data = result8, aes(x = iterazioni, y=learning_rate), color = "grey") 

Il valore del learning rate decresce all’aumentare delle iterazioni.

Cerchiamo un punto ottimale da dove partire per provare a cercare il punto di minimo globale (NB: questa tecnica non è sempre generalizzabile, va bene con questo tipo di funzione obiettivo): L’algoritmo calcola per ogni riga la funzione, dopo si prende le coordinate dove la funzione risulta più piccola, in modo tale da partire da un punto più basso e andare a cercare il minimo ottimale, anzichè scegliere un punto a caso come fatto in precedenza. Implementiamo solo l’algoritmo del gradient descent method e non quello stocastico perchè con questa quantità di dati è meno accurato.

findpoint <- function(coordinate){
  m <- as.matrix(coordinate)
  p <- c()
  pp <- c()
  pp2 <- c()
  x <- m[,1]
  y <- m[,2]
  for (i in 1:nrow(m)) {
    punto <- func_loss(c(x[i],y[i]))
    p[i] <- punto
    pp[i] <- x[i]
    pp2[i] <- y[i]
  }
  list(punto = p, coord = pp, coord2 = pp2)
}
c1 <- findpoint(coordinate)
c1 <- data.frame(minimo = c1$punto, x = c1$coord, y = c1$coord2)
subset(c1, minimo==min(minimo, na.rm=TRUE))

#Calcoliamo tre risultati diversi per scegliere il learning rate di partenza migliore 
res1 <- gdm(coordinate, c(308, 431), 0.1, 500, func_loss)
res2 <- gdm(coordinate, c(308, 431), 0.3, 500, func_loss)
res3 <- gdm(coordinate, c(308, 431), 0.5, 500, func_loss)
res1 <- data.frame(distanza = res1$minimo, iterazioni = res1$iterazioni, learning_rate =res1$learning_rate)
res2 <- data.frame(distanza = res2$minimo, iterazioni = res2$iterazioni, learning_rate =res2$learning_rate)
res3 <- data.frame(distanza = res3$minimo, iterazioni = res3$iterazioni, learning_rate =res3$learning_rate)
res1

ggplot() +
  geom_line(data = res1, aes(x = iterazioni, y=distanza, color = "red")) + 
  geom_line(data = res2, aes(x = iterazioni, y=distanza, color = "blue")) + 
  geom_line(data = res3, aes(x = iterazioni, y=distanza, color = "green")) + 
    labs(title="Andamento dell'algoritmo per ogni iterazione", 
        y = "distanza", 
        x = "numero di iterazioni") + 
    scale_color_manual("Learning rate \n di partenza", values = c("red", "blue", "green"), labels = c("0.1","0.3","0.5"))

Dai risultati del grafico emerge che partire da un learning rate di 0.3 si arriva più velocemente al punto di minimo ottimale. Partendo dalle coordinate (308, 431) sembrerebbe che l’algoritmo raggiunga un possibile punto di minimo globale, dove la la distanza è pari a: 2004.965.